2.6.8. Применение явной разностной схемы для решения
одномерного уравнения теплопроводности
Численные методы распадаются на два типа: один, явный, в котором
численное решение может быть выполнено шаг за шагом, исходя из
данного дифференциального уравнения и известных начальных и
граничных условий; второй, неявный, в котором неизвестные значения
связываются между собой системой линейных уравнений.
Рассмотрим применение метода конечных разностей для решения
следующей краевой задачи для уравнения теплопроводности. Пусть в
области
0 , 0x l t Q
требуется найти решение уравнения
2
2
2
,,
,
T x t T x t
a
t
x
(2.94)
удовлетворяющее начальному условию
,0 ( )T x F x
(2.95)
и граничным условиям
1
2
0,
( , ) ,
,T t f
T l t
t
ft
(2.96)
где
2
a
– некоторый постоянный коэффициент.
Если считать, что
2
a
есть коэффициент температуропроводности, то
решение этой задачи описывает нестационарное поле температуры
,T x t
в однородном стержне, теплоизолированном с боковой поверхности, при
заданном начальном распределении температуры
Fx
и температуре
1
ft
и
2
ft
соответственно на левом и правом концах стержня.
Если считать, что
2
a
есть коэффициент диффузии некоторой примеси
в сплошной среде, то решение этой задачи можно интерпретировать как
нестационарное поле концентрации примеси в трубке с известным
начальным распределением концентрации и заданными концентрациями
на левой и правой границах трубки.
Для решения этой задачи методом конечных разностей в области
изменения независимых переменных
x
и
t
построим прямоугольную
сетку с постоянными шагами
h
и
,
координаты узлов которой
определяются формулами (сетка и шаблон явной разностной схемы
представлены на рис. 2.10)
, 0,1, 2,..., ;
i
x ih i n
0,1, 2, ,..., .
j
t j j m
При этом узлы, лежащие на границе области
0 , 0 ,x l t
0, 0 ,x t Q
,0x l t Q
(обозначенные крестиками на рис. 2.10),
называются граничными узлами, остальные – внутренними.
Значения функции
,T x t
в узлах сетки будем обозначать символами
, .
j
i j i
T x Tt
Перейдем теперь к аппроксимации (во внутренних узлах
сетки) частных производных
,
ij
xt
T
t
,
2
2
,
ij
x t
T
x
в уравнении (2.94), а также
начальных и граничных условий в граничных узлах сетки конечно-
разностными соотношениями. Используя определение частной
производной по времени в точке
,
ij
x t
1
0
,
lim ,
ij
jj
ii
xt
T T T
t
x
0 = x
0
t
1
t
j
t
j + 1
t
m
= Q
x
1
x
i - 1
x
i
x
i +1
x
m
= l
t
h
Рис. 2.10
введем следующую аппроксимацию этой производной:
1
,
.
ij
jj
ii
xt
T T T
t
(2.97)
На основании определения частных производных по переменной
x
в
точке
,
ij
x t
соответственно справа
1
0
,
lim
ij
jj
ii
h
xt
T T T
x h
и слева
1
0
,
lim
i
j
jj
ii
h
xt
TT
h
T
x
для частной производной второго порядка по
x
получим
2
,,
2
0
,
lim
i
i j j
ij
t x t
h
x t
x
TT
xx
T
h
x
и введем для нее следующую аппроксимацию:
2
11
22
,
2
.
ij
j j j
i
x
ii
t
TTT
hx
T
(2.98)
Поскольку для аппроксимации частных производных использовались
значения искомой функции
,T x t
в четырех точках
1
, ,
ij
x t
,,
ij
x t
1
, ,
ij
x t
1
,,
ij
x t
то изображенная на рис. 2.10 жирными линиями
конфигурация узлов называется четырехточечным шаблоном.
Подставив аппроксимации (2.97) и (2.98) в уравнение (2.94) и
разрешив полученное уравнение относительно
1
,
j
i
T
получим разностное
уравнение
1
11
1 ,2
j j j j
i i i i
T T TT
(2.99)
2
2
1, 2,..., 1, 0,1, 2,..., 1, ,i n m
a
h
j
которое аппроксимирует исходное дифференциальное уравнение во
внутренних узлах сетки. Можно показать, что погрешность аппроксимации
имеет первый порядок по
и второй порядок по
,h
т. е.
2
.O h
Это
уравнение позволяет рассчитать значение функции
1j
i
T
на слое
1
1
j
jtt
по трем значениям
1
,
j
i
T
,
j
i
T
1
j
i
T
на предыдущем слое
.
j
t jt
В связи с этим реализуется следующий алгоритм послойного
расчета значений искомой функции в узлах сетки. Значения функции на
нулевом слое
0
0tt
вычисляются из начального условия (2.95):
0
, 0 ),( 0,1, 2,..., .
i i i
T T x F x i n
(2.100)
Во внутренних узлах первого слоя
1
tt
значения функции
1
, 1, 2,..., 1
i
T i n
рассчитываются по формуле (2.99). Недостающие
значения функции в граничных узлах первого и последующих слоев
получают из граничных условий (2.96):
1
0 1 1 1
1
1 2 1
0, ,
,.
n
T T t f t
T T l t f t
(2.101)
По найденному решению на первом слое аналогично находится
решение на втором и последующих слоях.
Совокупность разностных уравнений (2.99) и записанных в
дискретной форме начального (2.100) и граничных (2.101) условий
называется разностной схемой. Построенная разностная схема (2.99) –
(2.101) называется явной, поскольку в каждом уравнении (2.99) содержится
только одно значение
1j
i
T
неизвестной функции с последующего слоя
1
,
j
tt
которое определяется явным образом.
Достоинством этой схемы считается простота формул и алгоритма
вычислений. Недостатком является то, что ее можно использовать лишь
при выполнении условия
22
/ 0,5.ah
Это условие накладывает
сильное ограничение на величину шага по времени
,
т. е. вычисления
приходится вести с малым
,
что приводит к большим затратам
машинного времени. При использовании конечно-разностной схемы для
решения краевой задачи возникает вопрос об устойчивости такой схемы.
Под этим понимается следующее: конечно-разностная схема называется
устойчивой, если малые погрешности, допущенные в процессе решения,
затухают или по крайней мере остаются малыми при неограниченном
увеличении номера текущего слоя, в противном случае схема называется
неустойчивой. Явная схема является условно устойчивой.
Выведем далее формулы для расчета температуры в граничных узлах
при использовании граничных условий второго и третьего рода.
К граничным условиям второго рода приводят задачи, в которых
задана плотность теплового потока, протекающего через торцевые сечения
стержня. Пусть, например, через левое торцевое сечение стержня
поступает тепловой поток, плотность которого
1
0,qt t
есть заданная
функция времени. Приравнивая его к плотности теплового потока на левой
границе стержня, описываемой законом теплопроводности Фурье,
получаем
0
1
,
0, ,
t
q t t
T
x
где
– коэффициент теплопроводности материала стержня, т. е.
1
0,
.
t
T
x
t
(2.102)
Аппроксимируя производную в левой части конечно-разностным
соотношением, получаем
1
10
jj
j
t
T
h
T
или
0 1 1
.
jj
j
h
TtT
(2.103)
Аналогично, если через правое торцевое сечение стержня поступает
тепловой поток, плотность которого
2
,,q l t t
то
2
,
,.
lt
q l t t
T
x
(2.104)
Вновь аппроксимируя производную, имеем:
2
1
jj
j
nn
t
T
h
T
или
12
.
jj
n n j
h
TtT
(2.105)
Выражения вида (2.102) и (2.104), задающие
T
x
на границах стержня,
называются граничными условиями второго рода. Вытекающие из них
соотношения (2.103) и (2.105) позволяют рассчитать недостающие
значения функции в граничных узлах первого и последующих слоев.
К граничным условиям третьего рода приводят задачи, в которых на
торцевых поверхностях стержня имеет место теплообмен с окружающей
средой по закону Ньютона. В соответствии с этим законом плотность
теплового потока
,q
поступающего через торцевое сечение стержня,
пропорциональна разности известной температуры окружающей среды
c
T
и температуры
T
торцевого сечения стержня:
,
c
q TT
коэффициент пропорциональности
называется коэффициентом
теплообмена. Используя далее закон сохранения энергии и закон Фурье,
получим на левом торцевом сечении стержня соотношение
0,
0, .
l cl
t
T t T t
T
x
(2.106)
Аппроксимируя производную и переходя к принятым ранее
обозначениям температуры, имеем:
10
0
,
jj
j
l cl j
TT
TT
h
T
откуда
01
.
jj
l
cl j
ll
h
TTTt
hh
(2.107)
Аналогично, в случае теплообмена на правом торцевом сечении
стержня по закону Ньютона, получим соотношение:
,
,.
n cn
lt
T t T t
T
x
l
(2.108)
Снова аппроксимируя производную, имеем
1
,
jj
j
nn
n cn j n
TT
h
T t T
откуда
1
.
jj
n
n cn j n
nn
h
TtT T
hh
(2.109)
Выражения вида (2.106) и (2.108), задающие связь между
T
x
и
T
на
границах стержня, называются граничными условиями 3-го рода.
Вытекающие из них соотношения (2.107) и (2.109) позволяют рассчитать
недостающие значения функции в граничных узлах первого и
последующих слоев.
При формулировке краевых задач для уравнения теплопроводности на
каждой торцевой поверхности стержня можно ставить граничные условия
любого вида. Порядок точности разностного решения линейной задачи при
обеспечении условия устойчивости совпадает с порядком аппроксимации.
2.6.9. Применение неявной разностной схемы для решения
одномерного уравнения теплопроводности
Рассмотрим применение неявной разностной схемы для численного
решения уравнения теплопроводности (2.94) – (2.96):
2
2
2
,,
, 0 , 0 ,
T x t T x t
a x l t Q
t
x
,0 ( ),T x F x
1
2
0,
( , ) .
,T t f
T l t
t
ft
где
2
a
– некоторый постоянный коэффициент.
Для составления разностной схемы на сетке с постоянными шагами
h
и
по координатным осям
Ox
и
Ot
выберем изображенный на рис.
2.11 жирными линиями шаблон, состоящий из четырех точек
,
ij
xt
,
x
0 = x
0
t
1
t
j
t
j + 1
t
m
= Q
x
1
x
i - 1
x
i
x
i +1
x
m
= l
t
h
Рис. 2.11
11
,
ij
xt
,
1
,
ij
xt
,
11
,.
ij
xt
(Сетка и шаблон неявной разностной схемы
представлены на рис. 2.11.)
Входящие в уравнение (2.94) производные аппроксимируем
следующими конечно-разностными соотношениями:
1
1
,
,
ij
jj
ii
xt
T T T
t
(2.110)
1
2 1 1 1
11
22
,
2
.
ij
j j j
i i i
xt
T T T T
xh
(2.111)
Подставляя аппроксимации (2.110) и (2.111) в уравнение (2.94) и
перегруппировывая члены, получае:
1 1 1
11
1 2 ,
j j j j
i i i i
T T T T
(2.112)
1, 2,..., 1, 0,1, 2,..., 1,i n j m
где
2
2
a
h
.
Аппроксимация начального и граничных условий приводит к
соотношениям
0
,0 ( ), 0,1,..., ,
i i i
T T x F x i n
(2.113)
1
0 1 1
1
21
,
, 0,1, 2,..., 1.
j
j
j
nj
T f t
T f t j m
(2.114)
Построенная разностная схема (2.112) – (2.114) называется неявной,
поскольку в каждое уравнение получившейся системы линейных
уравнений (2.112) входят три значения неизвестной функции
111
11
,,
jjj
iii
TTT
с последующего слоя
1
1
.
i
j
tt
Можно показать, что
построенная неявная разностная схема имеет первый порядок
аппроксимации по
и второй по
.h
Реализуется следующий алгоритм послойного расчета значений
функции
,T x t
в узлах сетки. Значения функции
0
0,1,...,
i
T i n
на
нулевом слое
0
0tt
вычисляются из начального условия (2.113).
Расписав далее уравнение (2.112) для всех внутренних точек первого слоя
1
tt
и добавив соотношения (2.114) из граничных условий, получим
следующую систему линейных уравнений с известной правой частью для
определения неизвестных значений
1
0,1, 2,...,
i
T i n
на первом слое:
1
0 1 1
1 1 0
1 1 1
1
21
,
1 2 , 1, 2,..., 1,
.
i
i i i i
n
T f t
T T T T i n
T f t
По найденному решению на первом слое аналогично находятся
решения на втором и последующих слоях.
Преимуществом неявной схемы является то, что ее можно
использовать при любом значении
.
Неявная разностная схема считается
абсолютно устойчивой. Выбор шагов сетки определяется не
соображениями устойчивости, а требуемой точностью расчетов.
Недостатком неявной схемы является необходимость решать на
каждом временном слое систему линейных уравнений. Матрица этой
системы имеет «трехдиагональный» вид, поэтому для ее решения можно
использовать метод прогонки.
Выведем расчетные формулы метода прогонки. Идея метода прогонки
состоит в том, что путем последовательного исключения одного
неизвестного в каждом уравнении системы (2.112) ее можно привести к
эквивалентной системе с двумя неизвестными в каждом уравнении. Будем
искать решение системы (2.112) в виде
11
1 1 1
, 0,1, 2,..., 1,
jj
i i i i
T T i n
(2.115)
где
1i
и
1i
– неизвестные коэффициенты. Тогда
11
1
, 1, 2,..., .
jj
i i i i
T T i n
(2.116)
Подставим (2.116) в (2.112) получим
11
1
.
1 2 1 2
j
jj
ii
ii
ii
T
TT
(2.117)
Сравнивая формулы (2.115) и (2.117), выпишем рекуррентные
соотношения, связывающие коэффициенты
1i
и
1i
с
,:
ii
1
, 1, 2,..., 1,
12
i
i
in
(2.118)
1
, 1, 2,..., 1.
12
j
ii
i
i
T
in
(2.119)
Алгоритм решения системы линейных уравнений методом прогонки
состоит в следующем. Сравнивая выражение (2.115), записанное для
0i
11
0 1 1 1
jj
TT
, (2.120)
с граничным условием (2.114)
1
0 1 1
,
j
j
T f t
делаем вывод, что
1 1 1 1
0, .
j
ft
Далее по формулам (2.118), (2.119)
находим коэффициенты
1i
и
1i
1, 2,..., 1 .in
Этот процесс
называется прямой прогонкой. Зная прогоночные коэффициенты, решение
системы (2.112), (2.114) определяем по формуле (2.115) в обратном
порядке
1, 2,..., 2,1, 0 ,i n n
начиная с
1:in
11
1
,
jj
n n n n
TT
(2.121)
поскольку значение
1
21
j
nj
T f t
известно из граничного условия (2.114)
на правом конце отрезка. Этот процесс называется обратной прогонкой.
Выведем далее формулы для расчета
1
и
1
при использовании
граничных условий второго и третьего рода.
Пусть на левой границе стержня задано граничное условие второго
рода
1
0,
0,
t
t
T
qt
x
или в конечных разностях
11
0 1 1 1
.
jj
j
h
T T t
(2.122)
Сравнивая (2.122) с (2.120), делаем вывод, что
1
1,
1 1 1
.
j
h
t
Пусть и на правой границе стержня задано граничное условие второго
рода
2
,
,
lt
t
T
q l t
x
или в конечных разностях
11
1 2 1
.
jj
n n j
h
T T t
(2.123)
Решив систему двух уравнений (2.121) и (2.123), найдем, что
21
1
.
1
nj
j
n
n
h
t
T
Пусть на левой границе стержня задано граничное условие третьего
рода
л сл
0,
0,
t
T
T t T t
x
или в конечных разностях
11
л
01сл 1
лл
.
jj
j
h
T T T t
hh
(2.124)
Сравнивая (2.124) с (2.120), делаем вывод, что
1
л
,
h
л
1 сл 1
л
.
j
h
Tt
h
Пусть на правой границе стержня задано граничное условие третьего
рода
сп
,
,
п
lt
T
T t T l t
x
или в конечных разностях
11
1 сп 1
1.
jj
пп
n n j
hh
T T T t
(2.125)
Решив систему из двух уравнений (2.121) и (2.125), найдем, что
1
сп 1
1
j
пп
n j n n
hh
T T t
.
Метод прогонки можно применить, если знаменатель выражений
(2.118) и (2.119) не обращается в нуль. Это условие выполнено при
решении уравнения теплопроводности с граничными условиями всех трех
видов.
